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ABSTRACT 


In continuous-time optimal control problems, constraints must be satisfied as a set 
of logical conjunctions. In many practical space missions, however, the end-point 
functions may contain disjunctions. This thesis presents an approach for handling these 
end-point function disjunctions as part of a single continuous-time trajectory optimization 
problem. The approach embeds continuous representations of discrete logic operators as 
part of the problem formulation in order to model the disjunctions. The application of this 
new concept is first analyzed and illustrated for a canonical double integrator model as a 
proxy for practical space flight applications. It is then shown how the approach can be 
used to efficiently allow an algorithm to choose the minimum effort or minimum time 


attitude rotation for a rigid spacecraft amongst a set of possible terminal attitudes. 
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CHAPTER 1: 
Introduction 





1.1 Motivation 

Enabling spacecraft systems to more effectively plan and solve for minimum cost trajectory 
solutions can reduce the time and fuel required to maneuver between stations, increasing 
the number of potential collection events and fuel lifetime. Significant research has been 
done to develop autonomous trajectory and path planning solutions using Pontryagin’s 
Principle [1]-[4]. In these solutions path constraints and end-point conditions must be 
satisfied together as a conjunction. For example, it is common practice to define a desired 
end-point location in x, y and z coordinates AND a final velocity in x, y and z, AND path 


constraints. 


The purpose of this thesis work is ultimately to increase the mission value of on-orbit 
spacecraft. Many spacecraft missions are significantly limited by the ability to efficiently 
maneuver. Current solutions, to increase this efficiency, include a combination of simulation 
and human-in-the-loop commanding from the ground, often requiring multiple commands 
in sequence to achieve the desired end state. While significant work and research has 
been done to optimize and automate this process, expanding the current tools to include 
logical “OR” constraint conditions into the control and trajectory planning can significantly 
increase both spacecraft longevity and utility. A review of definitive texts on the subject 
including [1]-[4], and recent relevant work including [5]-[8] seem to validate that there 
is a gap in problem definition and solution development methodologies for problems with 


end-point disjunctions. 


1.2 Thesis Objective and Scope 

This thesis presents one methodology for implementing the Boolean exclusive or (XOR) 
operator into the end-point constraints for spacecraft slewing and pointing problems. This 
research first uses the two-dimensional double integrator model of a point mass as a sim- 


plified proxy for a spacecraft. A continuous representation of the XOR end-point condition 


is then defined. Once this new end-point definition is inserted into the two-dimensional 
optimal control problem (OCP), Pontryagin’s Minimization Principle is used to identify the 
necessary conditions for solution optimality, as well as any boundary condition needed to 
define an implementable boundary value problem (BVP). The problem is then solved and 


shown to meet the necessary conditions. 


This new concept for a disjunctive end-point constraint is then implemented using the 
quaternion representation of a rigid body rotation. This is done by first defining two math- 
ematical models (high and low fidelity) for spacecraft rotational dynamics and illustrating 
that these models well represent the actual systems they are modeling. Figure 1.1 shows a 
generic spacecraft in three-dimensional space, initially oriented so that the spacecraft body 
frame (red, green, blue) is aligned with the inertial frame (black). 

—o-X 


~o-Y 
—o-Z 


+2; +yi 


Figure 1.1: Initial spacecraft orientation in 3D space 


The XOR end-point constraint is then implemented to autonomously select between two 
different rotations, represented by two final quaternions, shown in Figure 1.2. Pontryagin’s 
Minimization Principle is again used to develop the necessary conditions, and the problem 


is solved and the solution shown to be optimal. 


Current solutions would require solving for each end-point, effectively enumerating through 
the disjunction, and then selecting the minimum cost option. This thesis work presents a 


novel and unique solution to this problem. 
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Figure 1.2: Spacecraft XOR reorientation problem 


1.3 Thesis Outline 


This thesis introduces the Lunar Reconnaissance Orbiter (LRO) mission as an example 
of a real world mission that can benefit from the introduction of disjuctive operators in 
optimal control. Chapter 3 introduces optimal control solution development using Pon- 
tryagin’s Minimization Principle. Chapter 4 details the development of the new method 
for implementing the XOR end-point condition. Chapter 5 then solves the OCP using the 
two-dimensional double integrator problem to demonstrate that this new technique is valid. 
This new methodology is then implemented for a spacecraft rotation using a quaternion 
representation. Chapter 6 constructs the dynamic model for this analysis, and then develops 
and solves the OCP using the more complex quaternion rotation model, demonstrating that 
this new technique can be expanded out to more complex space systems problems. Finally, 


Chapter 7 outlines the thesis conclusions and presents areas for future work. 
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CHAPTER 2: 
The Lunar Reconnaissance Orbiter Mission 





2.1 LRO Mission 


A timely example of how this research could impact real-world problems is the LRO 
mission. With NASA’s goal of returning humans to the moon by 2024, scientists need to 
learn as much as possible about the lunar surface, which makes the mission of the LRO 
vital. Improving the LRO’s collection efficiency would directly increase the amount of 
science information gathered during the resource constrained mission. The LRO is a robotic 
spacecraft that first entered lunar orbit on June 23, 2009, after making its way to space via 
an Alliance Atlas V rocket [9]. This started its Exploration Mission, focusing on furthering 


human existence in the solar system by mapping the moon. 


NASA has extended the LRO’s mission as it proved beyond capable of acquiring crucial 
data that addresses fundamental questions about lunar history and environment. Research 
questions included: 


Chronology/Bombardment, determining the bombardment history of the Moon 
from basin-forming events down to small, recent impacts; Crustal Evolution, 
investigating geological processes and their role in the evolution of the lu- 
nar crust and shallow lithosphere; Regolith Evolution, looking for evidence of 
processes that shaped the global lunar regolith in order to identify physical 
characteristics of the upper-most regolith layers; and Polar Volatiles, investigat- 
ing volatile sources, sinks and transfer mechanisms with emphasis on the lunar 


polar regions. [10] 


Using the above areas of research, scientists have used the LRO to draw important conclu- 
sions that will aid in future human habitation on the moon. In order to survive, colonists 
will need to have resources such as water readily available to them as cargo shuttles from 
Earth would be unable to keep up with the demand for necessary resources to sustain human 


life. Water maps and terrain maps created by the LRO will allow future settlers to pinpoint 
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hotspots for such resources and determine which regions would be of the greatest utility. 


Notable findings thus far include the coldest temperatures measured in the entire solar 
system in the polar shadowed regions of the moon, mammoth hydrogen deposits in regions 
where temperature supports the existence of water ice, a detailed index of surface craters as 
well as relatively new small impact areas, relatively new volcanic ranges, and unexpectedly 
large amounts of volatile gasses released near the lunar surface [9]. With the data collected 
via images taken by the LRO, NASA and supporting organizations have developed a 3-D 
map of the moon’s surface at 100-meter resolution with 98.2 percent coverage, excluding 
the polar shadow regions [11]. Shown in Figure 2.1 is an orthographic map of the Unified 
Geologic Map of the Moon that shows an updated topography of the lunar surface, developed 
by the USGS, NASA, and the Lunar Planetary Institute intended to be a reference for future 


moon missions [12]. 





Page = 


Figure 2.1: Geologic map of the moon. Source: [12] 


2.2 Instruments 

Such an in-depth mission requires highly-refined technology. The LRO has six highly so- 
phisticated instruments in its charge, as well as one technology demonstration, to include the 
Cosmic Ray Telescope for the Effects of Radiation (CRaTER), Diviner Lunar Radiometer 
Experiment (DLRE), Lyman Alpha Mapping Project (LAMP), Lunar Exploration Neutron 
Detector (LEND), Lunar Orbiter Laser Altimeter (LOLA), Lunar Reconnaissance Orbiter 
Camera (LROC), and Mini-RF Technology Demonstration [13]. 


CRaTER assesses the radiation environment on the moon to determine potential hazards 
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to human life and aid in developing technology to provide protection against such hazards. 
CRaTER is constructed of tissue-like plastic and measures the galactic and solar cosmic 
ray radiation to observe and analyze the effects that such radiation would have on human 
biology [14]. 


DLRE is the first instrument to measure the thermal emission from the moon’s surface. The 
instrument has measured the emissions during daily periods of shadow and illumination as 
well as during seasonal changes in shadow and illumination, which has revealed the moon’s 
extreme changes in temperature. Thermal imaging determines detailed temperature profiles 
on the moon, identifying cold traps and ice deposits that may endanger, or potentially 
benefit, human life. In addition to temperature hazards, DLRE also analyzes temperature 
differences to characterize topography, helping to identify potentially hazardous landing 
sites [15]. 


LAMP uses ultraviolet light to locate pockets of water ice that likely lie in the polar shadow 
regions. It uses UV light produced by stars as well as hydrogen atoms that exist sparingly 
in the solar system. In addition to finding these pockets of water ice, LAMP is also on a 
quest to pioneer Lyman-Alpha Vision Assistance (LAVA). LAVA uses natural starlight and 
ultraviolet skyglow to see in the dark, ultimately allowing astronauts to work in the darkness 
of lunar night and in the lunar poles [16]. 


Using thermal, epithermal and high-energy neutron sensors, LEND searches for hydrogen 
distributions one to two meters beneath the lunar surface in order to identify potential 
regions of water ice on the moon’s surface, particularly in the polar regions. The high 
resolution capability of LEND allows the detection of even the smallest distributions of 


hydrogen beneath the surface [17]. 


LOLA uses laser propagation and the subsequent range, scatter and return of the directed 
energy to determine the topographical profile of the moon. By identifying these elevation 
changes, scientists can make informed decisions regarding potential landing and exploration 
sites [18]. 


LROC consists of two narrow-angle cameras and one wide-angle camera that work together 
to “capture high resolution black and white images and moderate resolution multi-spectral 


images of the lunar surface” [10]. These cameras analyze data garnered via sunlight to 


provide information on the moon’s polar illumination conditions, thus helping to choose 
optimal landing sites. 


Finally, the Mini-RF Technology Demonstration consists of two radar instruments launched 
into space, the second of which resides on the LRO. The instrument is a radar suite used to 
locate water ice on the lunar surface. 


2.3. Mission Products 


Though many of these instruments appear to accomplish similar goals, the methods of 
data collection they use produce significantly different results and information. Figure 2.2 
provides images of the same lunar region, the South Lunar Pole, with very dissimilar results. 
Collectively, these instruments transmit 155 GB of data to Earth each day, all crucial to the 
future exploration and habitation of the moon. 
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Figure 2.2: LRO sensor modes. Source: [13] 


In addition to its primary mission mapping the moon, the LRO has been used to study the 
aftereffects of significant historical events, particularly lunar landings. Images of former 
Apollo landing sites have helped document the positions of many of the objects left on the 
moon from these missions, allowing scientists to confirm many assumptions made following 
the landings. The LRO first visited the site of the Apollo 11 landing, where humankind took 
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its first steps. The LRO captured images from just fifteen miles above the moon’s surface, 
and the images clearly depict objects from the mission, such as the lunar module and 
experimental instruments that the astronauts assembled, shown in Figure 2.3. Due to its 
high resolution capability, the LRO images even captured the tracks created on the lunar 


surface as the astronauts completed the first space walk [19]. 
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Figure 2.3: LRO image of Apollo 11 landing site. Source: [19] 


Following the visitation of the Apollo 11 site, the LRO re-visited the sites of Apollo 12, 
14, 15, and 16. During each of these visits, images from the LRO conveyed and confirmed 
details about the landing sites, including evidence of the first lunar rover’s exploration of 
the Mons Hadley massif near the Apollo 15 landing site [20] and the previously elusive 
Apollo 16 booster rocket impact site [21]. Not only do these revisits help showcase prior 
lunar activity, but they also allow scientists to compare data discovered during the Apollo 
Missions to data shown in current images. These comparisons help scientists make further 
hypotheses regarding the lunar environment. 


The Lunar Reconnaissance Orbiter plays a crucial role in future space exploration. Any 
errors in the LRO’s operation can result in, at best case, loss of important data or, at worst 
case, incorrect data that could one day risk lives. There are many natural phenomena that 
threaten the operability and capability of the LRO, so it is key that we maximize its utility 


while not providing additional threats from technologic error. 
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From the instruments on the LRO, scientists have produced the most sophisticated and 
accurate model ever of the lunar surface. The LRO needs to continue producing reliable 
images to support the endeavor of space exploration and needs to do so as efficiently as 
possible, to gather as much scientific data as possible in its finite lifetime. The satellite now 
uses a gyro-less control mode that prevents some collection targets from being accessed, 
if the star trackers are occulted [22]. Further development of controllers using algorithms 
with the non-conjunctive functionality presented in this work could automatically determine 
which of the possible collections are accessible without the need for enumerating through 
and trying them all. This would significantly increase the rate at which collection events 
were planned and increase the total number of science mission achieved. The LRO mission 
is an excellent example of the need for developing procedures for autonomously handling 
disjunctive constraint sets. Itis not the only space application that can benefit from a solution 
to this problem. Orbit transfers, landers, and mission planning are other examples where 


the ability to handle disjunctive constraints would be of benefit. 
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CHAPTER 3: 
Pontryagin's Minimization Principle in Optimal 
Control 





Optimal control problems are present in nearly every technical field. In these often complex 
problems, the goal is to develop a control solution that completes a task using minimum 
resource input. Intuitively, the optimal path is typically assumed to be a straight line. Though 
this can be the case for simple problems, “It has been shown for multi-body dynamic systems, 
like spacecraft antennas, that the optimal solution rarely follows a straight line because of the 
effects of nonlinearities and coupling between the various bodies” [6]. These problems are 
defined by their physical dynamics, written as a set of differential equations. The behavior 
of these systems is then further restricted by boundary conditions and constraints. Defined 
by Ross in his book, “a standard optimal control problem can be defined in terms of finding 
a dynamically feasible state-control function pair, x(-), u(-), that transfers the state system, 
x € R’*, froma given initial condition, x(t,) = x°, to a target condition e(x;, t+) = 0, while 
minimizing a given cost functional, J” [2]. Problem specifics vary widely with requirements 
and many other practical and engineering constraints, but irrespective of the details, the 
standard description of a continuous-time OCP is given as [2]: 


x € RX = (x1,...,xN,) ue R™ = (u,...,uy,) preamble 


Minimize: J[x(-),u(-), tf] = 


Osim E[Xo,Xpstote] +f’ Fx(t),u(t), dt } cost 

Control : Subjectto: x= f(x(t), u(t), t) } D 

Problem e! < e(x,,Xf,to tp) < e¥ rs) 
hh < h(x(t),u(t),t) < A" LH 


The cost functional J is the summation of an end-point cost functional EF and a running 
cost functional F. The inputs to this cost functional are the state variables (x), the control 
variables (u), time (t) and the final time value (f,). 
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The objective is to minimize J, subject to all of the linear and/or nonlinear constraints on 
the dynamics (D), end-point events (G), and path constraints (H/). For an optimal control 
solution, the constraints must be satisfied as a logical conjunction (see [23] for a discussion 
on the logical connectives). In other words, the logical expression D A & A H must be true 


for a feasible solution. 


In order to use common solution techniques such as collocation or shooting, the OCP may 
need to be further developed into a BVP. The BVP can be developed using Pontryagin’s 
Minimization Principle. While it does not solve the problem, Pontryagin’s Principle provides 
all of the necessary conditions for a candidate solution to be optimal. These are defined by 
Ross [2] in his Theorem 2.1: 


Theorem 2.1 (Pontryagin’s Principle) Given an optimal solution to Problem 
B, there exists an absolutely continuous covector function A(-) and a covector 


v that satisfy 
- the three Hamiltonian conditions; 
1. Hamiltonian minimization condition 
2. Hamiltonian value condition 
3. Hamiltonian evolution equation, 
- the adjoint equations and 


- the transversality condition. 


For a detailed discussion of the generalized steps in implementing Pontryagin’s Principle to 
develop BVPs and necessary conditions for optimality in OCPs as defined above, see [2], [6], 
[24], [25]. Two specific implementations of these generalized steps are presented here for a 
2D double integrator problem (Chapter 5) and for a 3D quaternion rotation (Chapter 6). 
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CHAPTER 4: 
Development of the Boolean XOR End Constraint 





In Chapter 3 the OCP was defined including a conjunction of problem dynamics (YD), end- 
point events (&), and path constraints (#1). But in many aerospace problems, the event set & 
may be formed by a logical expression involving non-conjuctive operations including nega- 
tions, disjunctions and/or exclusive-disjunctions across the end-point functions. Examples 
include selecting a feasible landing site from a set of alternatives, optimal reconfiguration 
of satellite formations, task selection and planning, or autonomous target engagement. An 
intuitive approach to solving such problems would be to enumerate across disjunctions. 
That is, first rewrite the original trajectory optimization problem as a collection of problems 
involving conjunctive constraints only. Then, solve the collection of problems to find the 


best solution. 


This chapter presents a new approach for handling disjunctive end-point conditions as part 
of a single continuous-time trajectory optimization problem without the need to resort to 
a potentially computationally intensive brute-force enumeration. The idea is to employ a 
continuous representation of the otherwise discrete logic network (see, for example, [26]) 
as a means to model disjunctive operators over the entire admissible problem state-space. 
This new concept is first illustrated using the canonical double integrator model as a suitable 


abstraction for an aerospace system. 


4.1 Continuous Representation of Discrete Logic Net- 
works 


An example logic network that can be found in many textbooks is shown in Figure 4.1. 
This logic network takes 3 inputs and implements the expression (x; V x2) A (2 V x3) by 
using the ‘AND’, ‘OR’, and ‘NOT’ operations which are normally instantiated in terms of 
electronic gates. Each of these elementary logical operations can also be represented by 
an appropriately defined continuous surface (see [27], [28]). These surfaces can then be 
combined in much the same way as the electronic gates of Figure 4.1 to perform the needed 
logic operation [26]. 
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Figure 4.1: Example of a typical logic network 


Amongst the various functional definitions that can be used to approximate discrete Boolean 
logic functions, radial basis functions (RBF) are perhaps one of the most simple, 

N 
(x) = >; aje Bllx-alll? 
i=l 

where WN is the number of RBF elements needed to generate a continuous surface that 
adequately approximates the discrete logic; c; locates the center of the surface peak 7, and a; 
and £ are appropriately chosen surface weights. Parameter 6 controls the slope of the RBF 
defined surface to transition between the binary logic states. The example logic network in 
Figure 4.1 can be implemented using the RBF parameters shown in Table 4.1, the output of 
which is shown plotted in Figure 4.2. 


Table 4.1: Radial basis function networks for basic logic functions. 




















operation | symbol RBF parameters 
AND =>- a= {+l}, c= {(1, 1)} 
B = variable 
D a= {+1,+1,+1}, c= {(1,0), (1,0), C1, 1)} 
OF B = variable 
XOR >>- a= {+1,+]}, e = {(1,0), (1,0)} 
B = variable 
S a= {+l}, c= {0} 
aoe B = variable 
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Figure 4.2: RBF representation of the network shown in Figure 4.1 for x3 = 1 and B = 12.5. 


4.2 Modeling a Two-Dimensional End-Point Disjunction 

This section details the combination of a set of two or more terminal conditions (x,y) 
into a single value for use as a continuous representation of the discrete exclusive or (XOR) 
logic. In practice, end-point constraints of this type can be implemented in coded algorithms 
using If = Then or Case logic. While this approach is easy to code, the results may be 


unpredictable because most optimal control solvers require a well-defined gradient. 


Another approach is to train a neural network to approximate the XOR function: 
e(xy,yf) = Net(xy, ys) =0 


To solve this problem, the neural network could be trained to solve an XOR problem (or 
any Boolean logic problem of any complexity), and the associated feed-forward equation 
of node weights, biases and activation functions can be extracted. This approach is im- 
plementable in a numerical solver because the gradients are the same as those needed for 
backpropogation during network training. Since backpropogation algorithms use differ- 
entiation, the extracted equations, while potentially complicated, can be differentiated to 


evaluate solution optimality. 


For the XOR case, one such neural network based construction is simply the summation 


of two RBFs as described in the previous section. Each RBF produces a Gaussian bump 
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with a peak value of one at its center and that tails sharply to zero as the distance from the 
center increases. Summing two of these functions, each centered at a valid point, results 
in a desirable surface. Mathematically, the end-point disjunction for a set of two possible 
end-points ((x¢,y¢) = (1,0) XOR (x, yf) = (0, 1)) can be modeled as: 


-[x)2+(y-1)?] -L-1)?+(y)?] 


f(x,y) =e 202 +e 202 


The surface has two RBF terms, one for each end-point. Inside each RBF is a summation 
of the squares of each term in the end-point. This being a two-dimensional problem, there 
are two terms in each RBF. This equation will expand for the more complicated quaternion 
(four-dimensional) case in Chapter 6. 


This particular form of the error surface places gaussian bumps at the points (0,1) and (1,0) 
—the desired end-points— with the surface gradient controlled by o. This allows for variation 
of the continuous representation to find the ideal compromise between approximating the 
discrete logic and the need for a smooth and continuous representation of the logic operator. 


Shown in Figure 4.3 is a representative surface using 7 = 0.3. 








Figure 4.3: Error surface for end-points (x,y) = (1,0) XOR (0, 1) 
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CHAPTER 5: 
Boolean XOR End Constraints in the 
Two-Dimensional Double Integrator Proxy Problem 





This chapter describes how to insert the logical XOR end-point constraint into a two- 
dimensional double integrator continuous-time optimal control problem formulation with 
consideration for the requirements of Pontryagin’s Minimization Principle. It then describes 
the process for implementing Pontryagin’s Minimization Principle to develop a BVP and 
solves the problem. 


5.1 Optimal Control Problem Definition and Analysis 

An idealized body subject to Newtonian physics is randomly placed in a two-dimensional 
plane. This body is modeled using the ubiquitous double integrator where * = uw; in this case 
¥ = ux, and ij = u,. Consider the problem,shown in Figure 5.1, of finding the optimal control 
that moves the body from its random initial rest position: x, y = (Xo, Yo) and Xo, Yo = (0,0) 
to rest at either valid end-point x,y = [(x¢,,yf,) ® (Xp.Yf))] and Xf, y¢ = (0,0) with 


minimum effort. 


Valid End 
Point 


N Random 
Starting 





Figure 5.1: Two-dimensional end-point XOR problem 


The desired 2-D double integrator system dynamics and required boundary conditions are 
then modeled as follows: 
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xX = Vy Kos Uo? S> ory”) 


System Vy = Uy Boundary VxooVy, = (0,0) 
Dynamics yo= Ww Conditions KPy¢ = PLO P2 
Vy = Uy Vxpo Vu = (0, 0) 


Optimal control of the double integrator is quite trivial (see [29]), but added complexity 
comes from the introduction of Boolean logic that defines a set of two acceptable boundary 
conditions. For initial conditions, the body is at rest at any random point in the 2-D space. 
The valid final position is defined as the exclusive or (XOR) of the two valid end-points. 
Using the RBF equation defined in the previous chapter and subtracting one to create a 
testable error equation results in the following end-point condition: 


=[(xp)?+(yp D7] -[0ep -1)?+(yp)?] 
Cxor(X fs, Yf) =e 2o2 +e 202 -—1 (5.1) 





As the search algorithm drives the point mass to either of the valid end-points, the value 


returned by Equation 5.1 would decrease to approximately zero. 


The dynamic optimization problem is now sufficiently described to allow for formulation 
and analysis using Pontryagin’s Principle. The dynamic optimization problem of interest is 


given in the standard form [2] as follows: 


se hele tis Aud _ 1 fit 
Minimize : J[x(-),u(-), tf] = 5 is (Crag Uy”) dt 
Subject to: ou = Dy, Vy 
Vx> Vy = Ux, Uy 
OCP: 
Vicor Vyor Vx ps Vyp = (0, 0, 0, 0) 
bask = (0, 5) 
—[oeg 2 +(yp 1)? ] —[org-D2 +p)? ] 
e 202 +e 20? =a 0 
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5.1.1 Developing the Necessary Conditions 

Along with developing the necessary conditions for optimality, solving this problem using 
shooting or collocation requires the formation of a BVP. Part of this process involves 
identifying additional boundary conditions as needed to properly form the BVP. These 
boundary condition are developed using Pontryagin’s Principle and the new XOR constraint. 
The Hamiltonian for the OCP is constructed [2]. 


2 4 2 
H(A,x,u,t) = F(x,u,t) + A’ f (x,u) = — + a FAYVy + AyxUy tAyVy + Avyly (5.2) 


Evaluating the adjoint equations (see [2]) by taking the partial derivative of the Hamiltonian 


(Eq. 5.2) with respect to each state variable describes the time variation of the costates. 


aH | 
ox Ay 0 A, is constant 
OH * 

OVx ze = Avx 2 Ax . Slope of Avx _ —A, 
oe Ay 0 A, is constant 
OH Avy Ay Slope of Avy = —Ay 
OVy 


Minimizing the Hamiltonian by taking the partial derivatives with respect to the controls 
defines the time-varying behavior of the controls, in this case equal to the opposite of their 
associated velocity costate. 

OH OH 


oe =Uy tay, = 0-0 uy, = Ay, Pm a ae Uy = —-Ay, 





Next, construct the Endpoint Lagrangian as in [2]. 
E(y,x(tp)) = E(x(tp)) + v2e(x(ty)) 
After inserting the terminal boundary conditions: 
Vp = Vy, = 0, tr=S, 


—[(xp)?+y¢-D?] —[(ep 1)? +(yp)?] 
€xor(Xf, Uf) =e 2o? +e 20? = 
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the Endpoint Lagrangian becomes: 


E(v, x,t) = Vo,(Vx~ — 0) + Vv, (Vy, — 0) + v1, (tf — 5). 


=[xep +p =] —[Org -D2 +p)? ] (5.3) 
t+ve |e 202 +e 202 =i 





The Hamiltonian Value Condition can be applied to determine the final Hamiltonian value 


by taking the partial derivative of the Endpoint Lagrangian (Eq. 5.3) with respect to tr. 


-0E 
H(t;) = Oty = S Vig 


The Hamiltonian should be equal to some unknown constant (v;, ) at ty. Taking the partial 
derivative of the minimized Hamiltonian with respect to time gives insight into how the 


Hamiltonian will evolve. 


dH OH _y 

dt Ot 
These last two conditions show that the Hamiltonian should be constant over the entire 
trajectory. 


Transversality analysis is used to determine the final values of the individual costates (A), 
which are equal to the partial derivative of the Endpoint Lagrangian with respect to the final 
states. 





A B 
JE =v -[(xp)?+(yp -D)?7] -[(xp -D?+(yp)?] 
A, (tf) = aes = —; (xf) e 22 +(x f = 1) e 202 ] 
Xf o : 
A B 
OE -y —[oxrp +p =] ~[xp D2 +(yp)7] 
Ay (tr) = ae 5 (yp—l)e 20? H(yp)e@ 20? | 
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By inspection, the position costates can be simplified. This is because for a given x, yr the 
complicated fractional terms are the same in each derivative, so they can be rewritten as 
arbitrary constants A and B as shown: 


Ax (tf) = —te (A * (xp) +B (xp - 1)) 


o2 


(5.4) 


The final position costates are equal to some unknown constant value and are coupled due 
to the disjunctive end-point error function. The final velocity costate values are equal to 


some arbitrary constants. 


OE 
Ay, (tp) = OvV;. = Vx ¢ 
Xf 


OE 
Ay, (tf) = dv, = Vyy 
if 


The coupled nature of the position costates can be leveraged to generate the needed event 
condition for formulation of the BVP. Rearranging the position costate equations (5.4) to 


isolate the y, coefficient terms: 


Ve. Ax (tf) 
o2 ((xp)(A +B) -B 


—Ve Ay (tp) 


a2 ((yp(A+B)-A 
After substituting, cross multiplying and rearranging: 
Ax (tf): [yp(A +B) — A] — y(t) - [xp(A + B) — B] =0 3:5) 


With the identification of an eighth event condition, the BVP is now sufficiently defined to 


allow coding and analysis via an algorithmic solver. 
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X,Y = Ux, Dy 





Vx> Vy = Ux, Uy 
Aehy = 0, 0 
Aves Ay = Ax, Ay 
XOR 
tos ty = (0, 5) 
BVP 
Xo> Yo = a5 y: 
~Leeg)?H(yp-DP1 Lp 1)? (yp)? ] 
e 202 +e 202 = 1 = 0 


(0, 0, 0,0) 
Ay(tp) « [xp(A + B) - Bl 


Vxo> Vyo> Vx > Vu 
Ax (te) » Lyp(A + B) — A] 


5.1.2 Simulation Results 

The solution to the XOR OCP — with an end-point disjunciton — is the minimum-effort 
state-control function pair tf — (x,u), which allows the point mass to be moved to the 
closer of the two possible end-points. In contrast to enumerating over the possibilities, the 
use of continuous representations of discrete logic operators allows the correct end-point 
to be chosen autonomously as part of the solution process. The inputs to the error surface 
function are the state variables xf and y,. 


In order to test the new concept of using a continuous representation of the discrete end-point 
disjunction, a series of dynamic optimization problems were solved. The first problem used 
a randomly generated starting point that was significantly closer to the (x, yr) = (0, 1) end- 
point. Such a point was chosen because it is known that the minimum effort solution would 
be to move the double integrator system to the point (0,1). The OCP with the end-point 
disjunction was solved using DIDO [30]. Twenty nodes were used because the solution, 
when propagated through the system dynamics, resulted in the correct end-point conditions 
(xf, yf) when rounded at the third decimal point. From the trajectory shown in Figure 5.2 


it is clear that the problem was properly solved with the selection of the closer end-point. 


Verification that the calculated solution is optimal, per the procedure outlined by Ross [2], is 


done by plotting the necessary conditions associated with Pontryagin’s Principle. Noting the 


pane 





[eos Trajectory 
xX Start 


© End 


> 0.5} 





0} O 


-0.5 - - 
-0.5 0 0.5 1 1.5 


[o, Yo] = (0.119, 0.498] —— > [x>, yr] = (0.000, 0.999] 


Figure 5.2: Example solution to the double integrator proxy problem showing that the correct 
solution was achieved. 


plot scale, Figure 5.3 shows that the Hamiltonian, as expected, is approximately constant. 
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Figure 5.3: Constant Hamiltonian evolution as required for optimality of the double integrator 
proxy problem. 


Shown in Figure 5.4, the minimization conditions, uy = —Ay, and uy = —Ay,, are met. 
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(a) Control trajectories (b) Velocity costates 
Figure 5.4: Control (wu) and velocity costate plots validating the u = —A, relationships for the 


double integrator proxy problem. 


Figure 5.5 shows that the 2, and A, costate trajectories are constant as expected and when 
compared with 5.4a they validate the Ay, + —A, and Ay, ~ —A, conditions. 
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Figure 5.5: Constant position costates equal to the slope of their control (uw) for the double 
integrator proxy problem. 


The position transversality condition (Eq. 5.5), using the algorithm values for final position 
and position costate, is approximately equal to zero (—2.3893E~°’), therefore all of the 


necessary conditions for optimality established using Pontryagin’s Principle [2] are met. 


To determine if a solution to the other valid end-point could be generated, a second problem 


was solved with a randomly generated starting point that was significantly closer to the 
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other end-point (xy, yf) = (1, 0). Examining Figure 5.6, it is again clear that the problem 
was properly solved with the selection of the closer end-point. 


— Trajectory 
x Start 


oO End 


> 0.5} 


-0.5 
-0.5 0 0.5 1 1.5 





x 


[9. Yo] = [0.856, 0.498] —— > ary, yy] = [1.000, 0.001] 


Figure 5.6: Example solution to the double integrator proxy problem showing that the correct 
solution was achieved. 


Figure 5.7 again validates that the Hamiltonian is approximately constant. 
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Figure 5.7: Constant Hamiltonian evolution as required for optimality of the double integrator 
proxy problem. 
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Shown in Figure 5.8, the minimization conditions, u, = —Ay, and uy = —Avyy, are again met. 
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Figure 5.8: Control (uw) and velocity costate plots validating the u = —, relationships for the 
double integrator proxy problem. 


Examining Figures 5.8a and 5.9, the A, and A, costate trajectories are again constant as 
expected, with Ay, ~ —A, and Ayy © —Ay. 
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Figure 5.9: Constant position costates equal to the slope of their control (uw) for the double 
integrator proxy problem. 


The position transversality condition (Eq. 5.5) was again calculated and is approximately 


equal to zero (7.7299E~®”). All of the necessary conditions for optimality established using 
Pontryagin’s Principle [2] are again met. 
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Defining which end-point is correct for given initial conditions is as simple as evaluating 
the Euclidean distance between the starting point and each end-point, and then selecting the 
shorter distance. In this case, that analysis divides the Cartesian plane along the diagonal 
line y = x, along which the distance between the two end-points in equal and propagating 


to either end-point results in an optimal solution. 


To illustrate how the correct solution can be found irrespective of the starting point, the 
simulation was repeated, generating random starting points throughout the input domain 
space. Figure 5.10 shows that for each initial position, the proposed OCP formulation that 
utilized a continuous surface to model the end-point disjunction made it possible to find the 
minimum effort end-point and generate a control trajectory solution that meets all necessary 


conditions for optimality. 
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Figure 5.10: Trajectory mapping illustrating the ability to select the correct end-point by modeling 
the end-point disjunction as a continuous function. 


pe 


5.2. Summary 

The simulation results presented in this chapter support the conclusion that this approach 
for implementing the exclusive disjunction (XOR) end-point constraint works as desired. In 
each case the coded algorithm was able to use the RBF defined error surface to drive the 
point mass to the optimal, minimum cost, end-point. The next logical step is to implement 


this end-point error surface in a more complex problem. 
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CHAPTER 6: 
Boolean XOR End Constraints in the Quaternion 
Reorientation Problem 





This chapter expands on the use of the end-point error surface for implementing functional 
logic in more complex optimal control problems. This is done by defining a dynamics model 
for spacecraft quaternion-based rotation, validating that model, and commanding rotations 
with the non-conjunctive end-point constraints. Control constraints limit four reaction wheel 
torques (AND), while an end-point event condition allows the final quaternion to be either of 
two valid options (XOR). The problem is then fully analyzed using Pontryagin’s Principle 
to identify the missing boundary condition and necessary conditions for optimality, and 
then coded in an algorithmic solver. It is then shown that the algorithm can independently 
identify the minimum time or minimum effort rotation and calculate an optimal control 


trajectory to rotate the model to that final quaternion. 


For this analysis the spacecraft body frame is initially aligned with the inertial frame. Control 
constraints on the torque applied to the reaction wheels effectively slow the spacecraft 
rotation and are intuitive to understand but difficult to visualize graphically. To begin to 
develop a valid solution for this scenario, a quaternion-based model for rigid-body spacecraft 
rotation is developed. This model is defined for spacecraft with four reaction wheels in the 
typical pyramid configuration aligned along the spacecraft yaw axis and uses differential 
equations for rotational dynamics and attitude kinematics. A series of tests are then presented 
to validate that the model well describes the rotational system. The entire constraint set is 
defined and Pontryagin’s Principle is then used to develop the necessary conditions for 
optimality. A coded algorithm is then used to determine optimal control trajectories for 
different combinations of rotations using the entire set constraints. The results are then 


analyzed to determine if the necessary conditions for optimality are met. 


6.1 Quaternion Model Development and Validation 
This section outlines the development of the spacecraft dynamics model using differential 


equations for attitude kinematics and rotational dynamics as defined in [8], [31]. This 
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is done using quaternions to define the spacecraft attitude with the following differential 


equation [8], [31], [32], such that the time rate of change of the quaternion is given as 
q = Q(w)q (6.1) 
Q(w) is the anti-symmetric matrix of spacecraft angular rates shown below. 


0 Wz Wy Wy 
1 |-w 0 Wy W 
Q(w) = 5 : 


Wy —a, O wr, 


—Wy Wy —Wz 0 


Euler’s equations for rotational motion are, [8] [31] 
I,-@ =—-Wsc X (Ise 5 ®h.,) = =: 


where I,,. is the spacecraft inertia matrix, w is the spacecraft angular acceleration and w 
is the spacecraft angular velocity. Additionally, h is the reaction wheel angular momentum 
and 7,3, is the reaction wheel control torque in the body-fixed coordinate frame, in a system 
with no external disturbance torques, henceforth referred to as spacecraft torque (T;-). All 


are expressed in vector form. Rearranging for spacecraft angular acceleration: 
. -1 
@=-I,.[w X (Iscw +h) - Tse] 


Since angular momentum is conserved for a zero net-bias system, the only torques on the 
body are those from the reaction wheels such that t,. = —Zt,, and the spacecraft angular 
acceleration simplifies to: 

Zt -w (6.2) 


@=-I 


Since h = Iw, the body-fixed reaction wheel angular momentum can be defined by multi- 
plying the reaction wheel angular rate (Q) and wheel inertia vector and then transforming 


that product into the body-fixed frame. 


5h, = Z(1,Q) 
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where Z is the matrix that rotates the reaction wheel values into the body-fixed coordinate 
system. Using the conservation of angular momentum, the derivative of the wheel angular 
momentum is equal to the sum of the torques on the body. 


ssa Mes) 
B w 
hy = 

dt 





= Zl yQy = —Ts = Zty 
Solving for the reaction wheel angular acceleration: 
Z1,Qy = Zt, 


0,215 7; (6.3) 


As discussed earlier four reaction wheel torques are modeled in a typical pyramid config- 
uration aligned with the yaw axis as in [8], [31] using the following Z matrix to apply a 


torque to the model, 


0.8192 0.8192 -0.8192 -0.8192 
Z= |0.4056 -0.4056 0.4056 —-0.4056 
0.4056 0.4056 0.4056 0.4056 


and the simplified inertia properties in the J matrix: 


Ixx 0O 0 
T5¢ = 0 lyy 0 kg : m- 
0 QO Izz 


The model is then defined by the state vector (Eqs. 6.1 - 6.3) and the control variable. 
ae 11 
x= 4 w Q| ER 
T 
u= [rr e Rt 


Assuming a zero-net bias system would allow for the use of the simpler form for the 


spacecraft angular acceleration (w), but for this relatively simple analysis the full version is 
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used. The model dynamics definition becomes: 


q = Q(w)q 
Problem 
@ =i} [w x (sew +h) — Tse] 
Dynamics 
QO =T;)u 


6.1.1 Model Validation 

To test the spacecraft rotational model, control torques (Figure 6.1) were defined and the 
system was propagated using Matlab’s ODE45 numerical integrator. Figure 6.2a shows, as 
expected, that any constant torque applied around an axis produced a linear angular rotation 


rate. 
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Figure 6.1: Quaternion rotation model - test control u = T 


Figures 6.2a and 6.2b show that a positive spacecraft angular rate causes the corresponding 


quaternion to move in the positive direction, and vice-versa. 


As always, changes to g4 should maintain the norm of the quaternion very close to one. 
Figure 6.3a validates that that is in fact the case as the norm error is very close to zero. 


Plotting the difference between the spacecraft and the reaction wheel angular momentum in 
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(b) Quaternion evolution 


Figure 6.2: Quaternion rotation model - spacecraft dynamic response 


each coordinate direction validates that angular momentum was conserved (Figure 6.3b). 
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(b) Conservation of angular momentum 


Figure 6.3: Quaternion rotation model - quaternion evolution and error norm 


The validation arguments support the conclusion that the model behaves as expected to 


simulate a rigid body rotating in space. Angular momentum bias is zero, so the simpler 


form for angular acceleration would be a good candidate for resource-intensive simulations. 


This model can now be used in solving constrained optimization problems. The next step is 


to identify the necessary conditions for optimality. 
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6.1.2 Optimal Control Problem Definition and Analysis 


The last requirement to finalize the quaternion rotation OCP is to define the XOR end-point 


constraint. 
qf (1) gp) 
pe OR) qf2(2) 
gf, (3) qf,(3) 
qf (4) qp,(4) 


Expanding the RBF surface discussed in Chapter 4 to a four-dimension end-point (four terms 
in each RBF), and again two valid end-points (two RBFs), the end-point error constraint 
(Eq. 6.4) becomes: 


-Uat ()-4 4, (D)? HaF (2)-a 4, DHF 3)-ap, 3))?+aF (4-45, (4))7] 





f e 20 a 
Ce -Uat (-ap (1)? +4F (2)-a 6 HGF 3)-ap, 8)? (aT 4-45, (4))? 1 (6.4) 
+e 202 — 1 





6.2 The Minimum Effort Problem 


Performing a full analysis via Pontryagin’s Principle is now possible using the quaternion 


XOR end-point constraint. A minimum effort cost functional is defined: [2] 


J[x(-),u(-),t¢] = E(x) + i: F(x, u)dt = sf (uy? aise +uz- +4”) dt 


to 
with the end-point constraints: 


Wy = wf =(0,0,0), 2% =Br = (0,0,0,0), qo = (0,0,0, 1) 


-LaF (ag, ()?+4F 2)-a5, (2)? CaF 3)-a 6, 3)? GT (4-44, (4)? I 
Cxor(q/) =e a eee 
-Lat ()-4p (1)? +(aF (2)-4 pg 2)? +aF (3)-ap, 3))?+(4F (4)-ap, (4)? 1 
+e 202 -1=0 








The Lagrangian of the Hamiltonian is constructed using the generic form (Eqn. 6.5). [2] 


A(p,A,x,u,t) = H(A,x,u,t)+uh(x,u,t) = F(x,u) +a" f(x,u) + yu h(x,u,t) (6.5) 
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= 2 2: 2 2 
H(p,A,x,u,t) = Ch eT 


+5 [Ag, (Gawx — 93W@y + 2Wz) + Ag, (G3Wx + G4Wy — Giz)... 
+4 [Ag (—qoWx + qiWy + G4Wz) + Agqy(—GiWx — G2Wy — 3Wz)]... 
Aw Ase [Z11U, + Zj2U2 + Z13U3 + Z14u4]... 


Awl 5em [Za1u1 + Z22U3 + Z23U3 + Zo4ua]... 
Aa Tec [Z3;U, + Z32U2 + Z33U3 + Z34U4]... 
+1Q, ees + Ao Ig, U2 + Aol os ,43 + Aas Tg},U4.-. 
FEU] + L2U2 + 3U3 + M4u4 
(6.6) 


Evaluating the adjoint conditions by taking the partial derivative of the Lagrangian of the 


Hamiltonian (eqn. 6.6) with respect to the states describes the time varying costate behavior. 


OH OH 1 








. 1 
8a1 = Aq, = zl Agnwz + Ag, Wy — Aqx] 0q3 Sag zag + Ag Wx — Aq] 
OH ; 1 OH ; 1 
Ian =-A,, = 5 Mae: — Ag,Wx — AgyWy] cen = —-Aq, = 5 aes t Ag Wy + Ags] 
OH : 1 
Jo =—-Ay, = 5 a4 + Aq93 — Ag392 - Aq, 91] 
OH ; 1 
dw, ten = 5 l-Aq 93 + gga + Aged ~ Aq? 
y 
OH ; 1 
A = Aw, = zag = Ap + G34 = Qq,93] 
Zz 
OH : OH : OH ‘ OH : 
dQ) 2 dQ, jos dQ, 


Often these equations are helpful in validating results, but not always. In this case the 
analysis has provided four easily testable, necessary conditions for the Ag behavior, namely 


that the values for Ag, are all constants. 


For the Lagrangian of the Hamiltonian to be minimized, both the stationarity and comple- 


mentarity conditions must be met: [2] 


OH 
Stationarity Condition : rr 0 (6.7) 
Uj 


35 


Control u=T a0 ay u; = —0.16 
Complementarity: pi-44=0 if -0.16<u;< 0.16 (6.8) 

Condition >0 if u;= 0.16 
These conditions are referred to as Karush-Kuhn—Tucker (KKT) conditions [2], [8], [31], 


and they generate two additional sets of conditions that can be easily tested. 


First, the stationarity condition can be graphed to show that the Lagrangian of the Hamilto- 
nian is constant over the trajectory. This is done by verifying that the derivative with respect 


to each control variable is mathematically zero over the trajectory. 


8 wy — AG TA, Zit — Auylod, Za — Au tse Za1 + Aa Tal + 11 = 0 

GH = yy — Ay TA, Zi2 — AwyloeyZ22 — Autres Z32 + Aaglal + lo = 0 a 
$F = ys — Ay TA, Z13 — Auylsd,Z03 — Aula, Z33 + Aaslgl + U3 =0 
BF = ug — Ay Tod Z14 — AuyloyZ24 — Aw Asi Z34 + Aagligh, + M4 = 0 





The complementarity condition shows a testable relationship between the control constraint 
and its 4 covector. If the control is at the upper limit, the covector should be positive, and if 
the control is at the lower limit, the covector should be negative. If the control is anywhere 
in between the limits, then the constraint is not active and the covector should be zero. 
Graphing the two together should show whether this relationship holds. 


Evaluating the partial derivative of the Endpoint Lagrangian (Eqn. 6.10) with respect to the 


initial and final times identifies the expected initial and final value of the Hamiltonian. 


E(y,x(tp)) = E(x(tp)) + v7 e(x(tp)) (6.10) 
E(v, x,t) = Vel€xor) + Vu, Wx + Vo,Wye + Va,Wzp 
oh fol f f y UE Ef (6.11) 
+V9,Q1 , + VQ, Qo, + ¥Q,Q3 + VQ, 24, + Vee (tf - t/) 
-0E 
On = H(t¢) =—Ver, 
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Taking the partial derivative of the Hamiltonian with respect to time: 


dH OH _ 
dt Ot 
These last two steps show that the Hamiltonian should be equal to some arbitrary constant 


over the entire trajectory. 


Transversality analysis is done to determine the final values of the individual costates, which 


are equal to the derivative of the Endpoint Lagrangian with respect to the final state. 


A 
Oo n‘£+£n+¢+¢§+¢+¢+ +++ +++ +++ + + 
-UaF (ag, (D)? +4F (2)-a4, (2)? a4 3)-a 6, 3)? HGF (4)-4 4, (4) | 
aglts)= Ze] PO 4H01¢ 
- OT 
-Uat ()-a 5 ())? HaF 2)-ap, DP HF 3)-ap, @))?HaT (4-45 (4)? 1 
tig! (i) -ap@le ce 








For brevity, each quaternion costate equation is expressed in terms of A and B, since for 


any combination of final quaternions, A and B are some constant value: 
Ve : ; ; . 
Ag lty) = ~-§ (Ala ~ an] + Bla ~ ap (0) (6.12) 


Rearranging (6.12) for the v. coefficient terms: 


—Ve Ag; (tf) 


o? — Alqf(i) — a7, 01 + Bla @ - ap] 





After substituting and rearranging, six new boundary conditions are developed: 


Cc 
_———SSSSSSSSSSS SS SS SSS Sh 


Agilt,)- (Ala! = an 1 + Bla! - ap) 


D 
 maaacnccncicncnnnncececcnnseen/~nseennencennnneeneneeenneneeeneneneneneeeen. 


= Ag (tp) (Ala!) - aq D1 + Bla!) - 4n)1) 


Vi, =[1,2,3,4], i #) 


(6.13) 
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Transversality analysis for the velocity costates is straightforward. 


OE 
Av(tf) = Aint = Vo 


OE 
Ag(tf) = 5Qf ~ 


The final velocity costate values are equal to some arbitrary constant which provides no 


VQ 


useful information. 


The BVP and the necessary conditions for solution optimality are now fully developed and 


the problem can be coded into an algorithmic solver. 


Minimize: J[x(-), u(), ty] 4 fF (DAL un?) dt 


Subject to: q = Q(w)q 
Oo = Tle x (sco +h) — tye] 
Q = I/u 
OCP: 
q° = [0, 0,0, 1]” 
exor(qh) = 0 
an (0, 0,0), (0,0, 0) 
Q,, Qf = (0, 0, 0, 0), (0, 0, 0, 0) 


6.2.1 Minimum Effort Results 

Once the minimum effort problem was fully defined and coded in an algorithmic solver, 
two different scenarios were presented and analyzed. Each minimum effort scenario was 
bounded to make the rotation in 300 seconds. Before using the algorithm to solve each 
XOR problem, the individual rotations were solved to identify baselines for cost, control 
trajectory and solution accuracy. This was done using the same dynamics model but with 


the standard 2, end-point event conditions: 
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q’.g! (0,0,0, 1), (q¢(1), a (2), 97 (3), 9 (4)) 
Wo, Wf = (0, 0, 0), (0, 0, 0) 
2,2; (0, 0,0, 0), (0, 0, 0, 0) 


Scenario 1: Angle Discrimination 

This scenario tested the algorithm’s ability to discriminate between two different magnitude 
rotations independent of any inertia tensor bias. It involves two end-point quaternions 
representing different magnitude rotations about different axes. To simplify the analysis, a 


spherically symmetric spacecraft body is assumed using a scalar inertia tensor matrix. 


sin(100/2) 0 
; 800 O O 
0 sin(40/2) > 
df= 0 ® , and Is-=| 0 800 O | kg-m 
0 O- 800 
cos(100/2) cos(40/2) 


The first quaternion represents a 100° rotation about the roll axis. Initial orientation and 
commanded rotation trajectories are shown in figures 6.4 and 6.5. Note that the control 
trajectory commands a pure rotation about the roll axis with a cost of 0.1628 kg?m7s. In 


this analysis, the kg*m*s have no actual context but are used as a relative measure. 


|—o =X 
|-o-Y 
|—O=Z 





Figure 6.4: Quaternion minimum effort angle discrimination, initial orientation for q, rotation 


20 


=—6 =X 
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Figure 6.5: Quaternion minimum effort angle discrimination, simple rotation to q, 


The second quaternion option represents a 40° rotation about the pitch axis. Initial orien- 
tation and commanded rotation trajectories are again shown in Figures 6.6 and 6.7. Again 
the optimal control trajectory commands a pure rotation about the pitch axis, this time at a 


much lower control cost of 0.1055 kg?m’s. 


—o-X 
—o-Y 


+Yyi 





Figure 6.6: Quaternion minimum effort angle discrimination, initial orientation for gz rotation 
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Figure 6.7: Quaternion minimum effort angle discrimination, simple rotation to q2 


Once these rotations were calculated, they were analyzed for optimality. Single end-point 
optimality analysis has been demonstrated in previous works (see [24]). That analysis is 
omitted here and instead a complete analysis of the XOR end-point case is presented. 


The minimum effort, angle discrimination, XOR problem was then solved. Figure 6.8 shows 
the initial spacecraft orientation with both possible end-point quaternions (g; and qz). Figure 
6.9 shows the resulting rotation. As desired, the algorithm correctly chose the minimum 
effort rotation (to gz) at a very similar control cost of 0.1054 kg?m7s. 

=-X 


—o-Y 
—O-Z 





+H 


Figure 6.8: Quaternion minimum effort angle discrimination, XOR rotation initial orientation 
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Figure 6.9: Quaternion minimum effort angle discrimination, XOR rotation 


It is interesting to note that there is a small difference in control cost between the two 


solutions. This is due to the difference in final end-point error and is discussed later. 


Optimality Validation 


Figure 6.10 shows that the adjoint necessary conditions are met as each 2,, is constant 


(A, = 0). The derivative of the Lagrangian of the Hamiltonian is approximately constant 


with respect to each control, validating the stationarity condition (Figure 6.10b). 
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Figure 6.10: Quaternion minimum effort angle discrimination, necessary conditions 
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The complementarity conditions (Figure 6.11), while met, are uninteresting as the constraint 


never activates. This is clear because the control costates (41;) never diverge from zero. 
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0.1} 0.5 0.1} 05 
hz hy 
: ee OT eee 
“ Mere 

-0.1} -0.5 -0.1} -0.5 
-0.2 . -1 0.2 — -1 

0 100 200 300 0 100 200 300 

Time (s) Time (s) 

oO Torque 3 Constraint 02 Torque 4 Constraint 

ao 0 SD OD OS. 0, OO ne 
0.1} 10.5 0.1} (0.5 

h, h 

: Opecse oe a eee 5a owen) - Ofeee se ews ae 0 6 o oogi) 

| bes bes 
-0.1} -0.5 -0.1} 1-0.5 

en ee a ie ar a apn ain ak a0 a 
-0.2 . . 1 -0.24 . . 1 

0 100 200 300 0 100 200 300 

Time (s) Time (s) 


Figure 6.11: Quaternion minimum effort angle discrimination, complementarity condition 


Figure 6.12 verifies that the Hamiltonian is constant throughout the trajectory. 
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Figure 6.12: Quaternion minimum effort angle discrimination, Hamiltonian 


Lastly, the transversality conditions were verified by calculating the four versions of equation 
6.13. In each case the difference was very close to zero, so all necessary conditions for this 


to be an optimal solution are met. 
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1.3059E~!2 
—1.2910E7!2 
3.1584E7!2 


Ag (t7) *C —Agg(t9) + D 
Ag (tf) uG — Aq, (tf) “D 
Ags (tf) -C —Agqy(tp) -D 


Scenario 2: Inertia Tensor Discrimination 

The second scenario tested the algorithm’s ability to discriminate between rotations of 
the same magnitude accounting for a non-symmetric inertia tensor. This case involves two 
possible end-point quaternions representing 40° rotations about different axes. A simplified, 


non-symmetric spacecraft body is assumed using a diagonal inertia tensor matrix. 


sin(40/2) 0 
712.55 O 0 
0 sin(40/2) > 
oF= ® and Is-=| 0 949.5 O | kg-m 
: 0 0 
0 0 871.5 
cos(40/2) cos(40/2) 


The first quaternion option now represents a 40° rotation about the roll axis. The commanded 

rotation trajectories are shown in figure 6.13. Note that the calculated control trajectory, as 

expected, commands a pure rotation about the roll axis with a cost of 0.0211 kg*m7s. 
[—o-X) 


|—@ = Y 
j—O =Z 





yf SS. OS 


y x 


Figure 6.13: Quaternion minimum effort inertia discrimination, simple rotation to q, 


As in Scenario 1, the second quaternion option represents a 40° rotation about the pitch 
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axis. The commanded rotation trajectories are again shown in figure 6.14. Again the optimal 
control trajectory commands a pure rotation about the pitch axis, this time at a much higher 
control cost of 0.1496 kg?m?s. 


-—o-X 


[ro =2) 


+Yyi 





Figure 6.14: Quaternion minimum effort inertia discrimination, simple rotation to q2 


The minimum effort, inertia tensor discrimination, XOR problem was then solved. The 
algorithm correctly chose the minimum effort rotation again at a very similar cost to the 
single end-point rotation problem of 0.0205 kg*m?s. 

=x] 


j= ~Y 
0-7 


+41 





Figure 6.15: Quaternion minimum effort inertia discrimination, XOR rotation initial orientation 
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Figure 6.16: Quaternion minimum effort inertia tensor discrimination, XOR rotation 


Optimality Validation 
Figures 6.17 and 6.18 show that the adjoint necessary conditions, stationarity, and constancy 
of the Hamiltonian are all met. The torque constraints remain inactive so the complemen- 


tarity conditions are again uninteresting. 
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Figure 6.17: Quaternion minimum effort inertia tensor discrimination, necessary Conditions 
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Figure 6.18: Quaternion minimum effort inertia tensor discrimination, Hamiltonian 


The transversality conditions were again verified to be very close to zero. 


Ag (tp)-C—-Ag(ts)-D = 5.3693E-% 
Agltf)-C —Ag (ts) D = 2.9692E~% 
Ag (tp) C —Ag(tf)-D = —4.9017E~% 


Again all necessary conditions for this to be an optimal solution are met. 


6.3 The Minimum Time Problem 


This section presents the analysis of the same XOR problem, now for a minimum time cost 
functional [2] 


J[x(-),u(-),t¢] = E(xyz) + o F(x,u)dt = ty (6.14) 


with the same end-point constraints: 


Wo = wy = (0,0, 0), OQ, = QA = (0,0, 0, 0), do = (0,0, 0, 1) 


-(aF (ag, ())? af (2)-a5, (2)* CaF 3)-a 4, 3) GT (4)-a 4, (4)? ] 








Cxor(q/) =e 202 ee 
-Lat ()-4ap (1)? +(aF (2)-4 pg (2)°+aF (3)-ap, 3))? CaF (4)-ap, (4)? ] 
+e 202 =-1=0 
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Using equation 6.5, noting the loss of the u* terms, the Lagrangian of the Hamiltonian is 


now linear in control [2], thus we expect a band-bang type control. 


A(p, A,x,u,t) = 51g, (q4wx — G3Wy + q2Wz) + Ag (Q3wx + G4Wy — qiWz) |... 
+5 [Ag3(—42Wx + F1Wy + Gawd) +Agy(—G1Wx ~ G2Wy — G302)]... 
=Ais Teen [Z11u1 + Zi2U2 + Z13uU3 + Z14u4)... 
Aw Ten [Zo1U, + Z22U3 + Z23U3 + Zo4ua)... 
Aw Tyey,(Z3iu + Z32u2 + Z33U3 + Z34ug)... 


-1 -1 -1 -1 
Aa lin, 41 + Aaa ly, 42 + Aas lin;343 or Rolla 


FU] + L2U2 + (3U3 + L4u4 
(6.15) 


Evaluating the adjoint conditions results in the same testable conditions as before. Likewise, 
the same KKT conditions (stationarity Eq. 6.7 and complementarity 6.8) must be met [2]. 


The stationarity condition equations no longer contain the control term (u). 


gH 7 Aa Lsey)ZA1 ae ee Z21—A Tye 231 + Aa Ip), + uy =0 


L Wy" $c22 Wz S033 
GH = — Aer Lec, Z12 — Ay gy 222 — Aig] seyg 232+ AQgT gt, + U2 = 0 6.16) 
GH = — Aer Lsey, Z13 — Atay] gey223 — Ate] sing 233 + ADs lat, + U3 = 0 
ou = Ay Tip, 214 — Aeyl sey 224 — Aw AD 5ey32Z34 + Agalg,, + M4 = 0 


The complementarity condition plots should still show the relationship between the control 
constraint and its 44 covector. Unlike in the previous case, for the minimum time problem 


the constraints should become active. 


The Endpoint Lagrangian now includes a ty term. 


E(y,x(tp)) = E(x(tp)) + v7 e(x(ty)) (6.17) 


E(vx7,t7) = tp + Ve(€xor) FV aga, + VoWyp + Vin, Wz (6.18) 
+V¥9, 24 , + VQ, 9229 + ¥9,3 , a VQ4224 


WE dH oH 
oe eye Ae did ee 
(7) UES Sa Bs 


0 
Ot f 
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These last two equations now show that the Hamiltonian should be equal to negative one 
and constant over the entire trajectory. Transversality analysis is unchanged with the same 


quaternion costate behaviors (6.13). 


The OCP and the necessary conditions for solution optimality are again fully developed and 


the problem can be coded into an algorithmic solver. 


Minimize: J[x(-),u(-),t¢] = tf 
Subject to: q = Q(w)q 
@ = -Iyt[o x (yew +h) - tc] 
OCP: . : Tu 
q° = [0,0, 0, 1]? 
exor(q’!) = 0 
Wo, WF = (0, 0,0), (0, 0, 0) 
Q,, Qf = (0, 0, 0,0), (0, 0, 0, 0) 


6.3.1 Minimum Time Results 
This scenario tested the algorithm’s ability to solve for a minimum time rotation by choosing 
between different magnitude rotations. This test again involves two end-point quaternions 


representing different magnitude rotations about different axes. 


sin(120/2) 0 He G 0 
0 0 2 
df= ® and I;-=| 0 949.5 O | kg-m 
0 sin(30/2) 
0 0 871.5 
cos(120/2) cos(30/2) 


The first quaternion represents a 120° rotation about the roll axis. Initial orientation and 
commanded rotations are shown in figures 6.19 and 6.20. Note that the control trajectory 


commands a pure rotation about the roll axis with a rotation time of 107.40 seconds. 
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Figure 6.19: Quaternion minimum time, initial orientation for q, rotation 
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Figure 6.20: Quaternion minimum time, simple rotation to q1 


The second quaternion represents a 30° rotation about the yaw axis. Initial orientation and 
commanded rotation trajectories are again shown in figures 6.21 and 6.22. Again the optimal 
control trajectory commands a pure rotation about the pitch axis with a rotation time of 
84.37 seconds. 


As in the previous, minimum effort case, the baseline control trajectories and rotations were 
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Figure 6.21: Quaternion minimum time, initial orientation for gz rotation 
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Figure 6.22: Quaternion minimum time, simple rotation to q2 


calculated and then verified to meet all necessary conditions for optimality. That analysis is 


again omitted in favor of a complete optimality analysis of the XOR problem. 


The XOR problem was then solved. Figure 6.23 shows the initial spacecraft orientation with 
both possible end-point quaternions (qg; and q2). Figure 6.24 shows the resulting rotation. 
As desired, the algorithm correctly chose the minimum time rotation (to g2) at a very similar 


total rotation time of 84.30 seconds. 
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Figure 6.23: Quaternion minimum time, initial orientation for the XOR rotation using non- 
symmetric inertia tensor 
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Figure 6.24: Quaternion minimum time, XOR rotation using non-symmetric inertia tensor 


Optimality Validation 
Figure 6.25 shows that the adjoint necessary conditions are met (each 2g is approximately 
constant, Ag = 0), and that the derivative of the Lagrangian of the Hamiltonian is approxi- 


mately constant with respect to each control, validating the stationarity condition. 


This complementarity condition says that when path constraint (7) is at the upper constraint, 


the covector must be positive or zero. Inversely, if f is at the lower constraint the covector 
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Figure 6.25: Quaternion minimum time, necessary conditions 


must be negative or zero. If fis in between, the constraint is inactive and the covector should 
be zero. In each constraint case, shown in Figure 6.26, complementarity is met as expected 


for the minimum time bang-bang type control. 
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Figure 6.26: Quaternion minimum time, complementarity condition validation 
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Figure 6.27 verifies that the Hamiltonian is constant and approximately negative one 
throughout the trajectory, with the slight bump in the middle as commonly observed in 


numerical solutions to minimum time OCP. 





Hamiltonian 








0 20 40 60 80 100 
Time (s) 


Figure 6.27: Quaternion minimum time, Hamiltonian 


Lastly, the transversality conditions were verified to again be close to zero. 


Ag (tf) °C —Ag (tf) D = -2.8820E-” 
Agltf)-C-Ag(ts):D = 0.0541 
Ag, (tf) -C = Aas te) -D = 0.0361 


These transversality condition values are much higher than in the previous minimum effort 
problem. While larger, the final error is still small, and comparing the XOR trajectory to 
the simple rotation for the same end-point it is clear that the XOR solution is the same and 


therefore is evaluated as meeting all necessary conditions to be an optimal solution. 


6.4 Summary 

The simulation results presented in this chapter again support the conclusion that the newly 
developed approach for implementing the exclusive disjunction (XOR) end-point constraint 
can be applied for practical space problems —in this case an optimal spacecraft reorientation. 
In each case the coded algorithm was able to use the developed continuous error surface to 
choose the minimum cost rotation and develop a control solution to rotate the spacecraft to 


the correct final quaternion state, while meeting all the necessary conditions for optimality. 
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The goal of this Chapter was to implement a method for choosing between two reorientations. 
It did not specifically address how well the algorithm reached the final state. It is however 
important to note, as highlighted earlier, that there were very small differences between 
the calculated trajectories, and final cost, between the XOR rotations and the single end- 
point rotations. This was due to small differences in the final end-point error. The more 
complicated end-point definition reduced the final accuracy very slightly. 


For the minimum time rotation the final quaternion target was 


0 
ie: 
4” ~ 10.258819 

0.965925 


while the single end-point rotation ended at 


9.0058E~7 

_ |-8.8416E~"° 
ISR) 0.2588 
0.9659 


and the XOR rotation ended at 


T1059E-8 
—1.4509E77 
0.2584 
0.9660 


dXOR = 


These very small differences in the final end-point were, in the most significant case, at the 
fourth decimal. These differences resulted in correspondingly small differences in trajectory 
and final cost. While the difference is small, it is not zero. It is important to consider what 
degree of final position accuracy is required for a given implementation. For a satellite 


slewing maneuver, this solution is sufficient to be implemented in a planning architecture 
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to determine the minimum cost rotation. If it were to be implemented in an actual slewing 
controller, the XOR solution is accurate enough to rotate a spacecraft to the point that a 
feedback controller could drive the residual error to any needed requirement. 
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CHAPTER 7: 
Conclusions and Future Work 





7.1 Conclusions 

This thesis outlined the need to include non-conjunctive constraints as part of continuous- 
time optimal control problems. Expanding the available logic connectives for OCPs beyond 
conjunctions, can enable spacecraft like the Lunar Reconnaissance Orbiter to more effec- 
tively plan and solve for minimum cost trajectory solutions. This would reduce the time and 
fuel required to maneuver between stations and increase the number of potential collection 
events, fuel lifetime, and total mission value. While significant research has been done 
to develop autonomous trajectory and path-planning solutions, this thesis expands on that 
body of work, detailing one method for implementing logical XOR end-point constraint 


conditions into the control and trajectory planning. 


The two-dimensional double integrator model of a point mass is a reasonable, albeit sim- 
plified, proxy for a spacecraft. This initial piece of the thesis work effectively demonstrated 
that, in this idealized application, this new method can be used to effectively implement the 
XOR end-point constraint. Expanding the implementation of this new method to a more 
complex spacecraft reorientation model was an integral step in the development of this 
new optimal control tool. This work was able to effectively demonstrate that a trajectory 
optimization solver could autonomously determine which final quaternion represented the 


minimum cost rotation. 


7.2 Future Work 


This work should continue to evolve. Implementation in increasingly complex systems and 
eventually hardware is a logical and achievable goal. One interesting topic to explore going 
forward is the possibility of inserting parameters into the problem formulations (see [30]). 
It is likely that the use of end-point parameters could reduce the final error when using more 


complex end-point error functions. 


Another avenue for future work to continue exploring the question of non-conjunctive 
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constraints is the use of non-conjunctive logic in trajectory path constraints. There is a 
subtle but important difference between end-point and path constraints that alters the way in 
which they are implemented. That is that end-point constraints are valid only for an instant 


in the trajectory, while path constraints must be true throughout. 


Additionally, exploring logical operations other than OR/XOR could prove valuable. Devel- 
oping a method for implementing the Boolean NAND (not AND) block would be particularly 
valuable. This is because in Boolean logic the NAND block has the property of functional 
completeness, meaning that it is possible to develop any truth table using only the NAND 
blocks. If a method for implementing the NAND block is developed, then it is possible 
that any Boolean truth table function could be implemented into continuous-time optimal 
control problems in the manner detailed in this thesis. 


58 





List of References 





[1] L. S. Pontryagin, Mathematical Theory of Optimal Processes. Routledge, 2018. 


[2] I. M. Ross, A Primer on Pontryagin’s Principle in Optimal Control, 2nd ed. Colle- 
giate Publishers, 2015. 


[3] D. E. Kirk, Optimal Control Theory: An Introduction. Courier Corporation, 2004. 


[4] A. E. Bryson, Applied Optimal Control: Optimization, Estimation and Control. 
Routledge, 2018. 


[5] P. C. Calhoun and J. C. Garrick, “Observing mode attitude controller for the lunar 
reconnaissance orbiter,” in Proceedings of the 20th International Symposium on 
Space Flight Dynamics, 2007. 


[6] R. Arledge, “Implementation of optimal controls using conventional control sys- 
tems,” Master’s thesis, Dept. of Mech. and Aerosp. Eng., Naval Postgraduate School, 
Monterey, CA, 2014. 


[7] J. Kaufman, “Automated maneuver design and checkout for the Lunar Reconnais- 
sance Orbiter,’ Master’s thesis, Dept. of Mech. and Aerosp. Eng., Naval Postgradu- 
ate School, Monterey, CA, 2014. 


[8] T. Lippman, “Enhancing the science collection capability of NASA’s Lunar Recon- 
naissance Orbiter (LRO),” Master’s thesis, Dept. of Mech. and Aerosp. Eng., Naval 
Postgraduate School, Monterey, CA, 2017. 


[9] The LRO Mission. NASA. [Online]. Available: https://lunar.gsfc.nasa.gov/about. 
html. Accessed Apr. 20, 2020. 


[10] LROC. [Online]. Available: https://www.lroc.asu.edu. Accessed Apr. 20, 2020. 


[11] NASA Probe Beams Home Best Moon Map Ever. (2011, Nov. 18). SPACE.com. 
[Online]. Available: https://www.space.com/13666-moon-map-lunar- 
reconnaissance-orbiter.html. Accessed Apr. 20, 2020. 


[12] USGS Releases First-Ever Comprehensive Geologic Map of the Moon. (2020). 
USGS. [Online]. Available: https://www.usgs.gov/news/usgs-releases- first-ever- 
comprehensive-geologic-map-moon. Accessed Aug. 19, 2020. 


[13] The LRO Instruments. NASA. [Online]. Available: https://lunar.gsfc.nasa.gov/ 
instruments.html. Accessed Apr. 20, 2020. 


59 


[14] 


[15] 


[16] 


[17] 


[18] 


[19] 


[20] 


[21] 


[22] 


[23] 
[24] 


[25] 


Cosmic Ray Telescope for the Effects of Radiation. (2015). University of New 
Hampshire. [Online]. Available: http://crater.sr.unh.edu/. Accessed Apr. 20, 2020. 


D. Paige and B. Greenhagen, “Diviner lunar radiometer experiment extended mis- 
sion results: Thermal, thermophysical and compositional properties,” in Lunar and 
Planetary Science Conference, 2013, vol. 44, p. 2492. 


The Lyman Alpha Mapping Project: Seeing in the Dark. Southwest Research In- 
stitute (SwRI). [Online]. Available: https://www.boulder.swri.edu/lamp/. Accessed 
Apr. 20, 2020. 


I. Mitrofanov, A. Bartels, Y. Bobrovnitsky, W. Boynton, G. Chin, H. Enos, L. Evans, 
S. Floyd, J. Garvin, D. Golovin et al., “Lunar Exploration Neutron Detector for the 
NASA Lunar Reconnaissance Orbiter,’ Space science reviews, vol. 150, no. 1-4, pp. 
183-207, 2010. 


Lunar Orbiter Laser Altimeter. NASA Goddard. [Online]. Available: https://lola. 
gsfc.nasa.gov/. Accessed Apr. 20, 2020. 


Apollo 11: ‘A Stark Beauty All Its Own’. NASA. [Online]. Available: https://www. 
nasa.gov/mission_pages/LRO/news/apollo-11.html. Accessed Apr. 20, 2020. 


Apollo 15: Interplanetary Mountaineers. NASA. [Online]. Available: https://www. 
nasa. gov/content/goddard/apollo- 15-original-interplanetary-mountaineers. Accessed 
Apr. 20, 2020. 


LRO Finds Apollo 16 Booster Rocket Impact Site. (2019, Apr. 3). NASA. [Online]. 
Available: https://www.nasa.gov/image-feature/goddard/lro-finds- apollo- 16-booster- 
rocket-impact-site. Accessed Apr. 20, 2020. 


J. Halverson, P. Calhoun, O. Hsu, J.-E. Dongmo, R. Besser, B. Ellis, R. De Hart, 

Y. Tedla, S. Rosney, and S. T. Snell, “Testing of the lunar reconnaissance orbiter at- 
titude control system re-design without a gyro,” in 42nd Annual AAS Guidance and 
Control Conference. Breckenridge, CO: AAS 19-108, 2019. 


E. Mendelson, /ntroduction to Mathematical Logic. CRC Press, 2009. 


A. Fleming, P. Sekhavat, and I. M. Ross, “Minimum-time reorientation of a rigid 
body,” Journal of Guidance, Control, and Dynamics, vol. 33, no. 1, pp. 160-170, 
2010. 


A. Fleming, “Real-time optimal slew maneuver design and control,” Master’s thesis, 
Dept. of Mech. and Astro. Eng., Naval Postgraduate School, Monterey, CA, 2004. 


60 


[26] S. Shi, H. Chen, M. Zhang, and Y. Zhang, “Neural logic networks,” arXiv preprint 
arXiv: 1910.08629, 2019. 


[27] S. Haykin, Neural Networks: A Comprehensive Foundation, 2nd ed. Prentice-Hall, 
Inc., 2000. 


[28] C. M. Bishop et al., Neural Networks for Pattern Recognition. Oxford University 
Press, 1995. 


[29] M. Athans and P. Falb, Optimal Control an Introduction to the Theory and Its Appli- 
cations. McGraw-Hill, 1966. 


[30] I. M. Ross, “Enhancements to the DIDO optimal control toolbox,” arXiv preprint 
arXiv: 2004.13112, 2020. 


[31] T. Lippman, J. Kaufman, and M. Karpenko, “Autonomous planning of constrained 
spacecraft reorientation maneuvers,” in AAS/AIAA Astrodynamics Specialist Confer- 
ence. Stevenson, WA: AAS 17-676, 2017. 


[32] M. J. Sidi, Spacecraft Dynamics and Control: A Practical Engineering Approach. 
Cambridge University Press, 1997. 


61 


THIS PAGE INTENTIONALLY LEFT BLANK 


62 





Initial Distribution List 





1. Defense Technical Information Center 
Ft. Belvoir, Virginia 


2. Dudley Knox Library 


Naval Postgraduate School 
Monterey, California 


63 


